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SUMMARY 


A "zonal," or "patched," grid approach is one in which the flow region of inter- 
est is divided into subregions which are then discretized independently, using exist- 
ing grid generators. The equations of motion are integrated in each subregion in 
conjunction with zonal-boundary schemes which allow proper information transfer across 
interfaces that separate subregions. The zonal approach greatly simplifies the 
treatment of complex geometries and also the addition of grid points to selected 
regions of the flow. In this study a conservative, zonal-boundary condition that 
could be used with explicit schemes has been extended so that it can be used with 
existing second-order-accurate implicit integration schemes such as the Beam-Warming 
and Osher schemes. In the test case considered, the implicit schemes increased the 
rate of convergence considerably (by a factor of about 30 over that of the explicit 
scheme) ♦ Results demonstrating the time-accuracy of the zonal scheme and the feasi- 
bility of performing calculations on zones that move relative to each other are also 
presented. 

INTRODUCTION 


A major problem in computational fluid dynamics is the generation of grids for 
realistic, three-dimensional bodies such as complete aircraft configurations. The 
generation of a single grid that discretizes the entire flow region associated with a 
complex configuration is extremely difficult and, for some problems, impractical. 
Often, an attempt to generate a single grid for a complicated flow region results in 
highly skewed grids which in turn result in inaccurate calculations; also, such an 
attempt usually requires an inordinate amount of the aerodynamicist ! s time to "tune" 
existing grid-generation schemes to yield acceptable grids. A second factor that 
contributes to the complexity of grid generation is the necessity to cluster grid 
points in regions where the dependent variables and their gradients change rapidly 
(selective grid refinement) . 

The problems mentioned above can be overcome to a limited extent by developing 
very sophisticated grid-generation schemes. However, an alternative approach is to 
divide the given region into simpler subregions such that each subregion (or zone) is 
a geometrically simple figure; for example, every two-dimensional region can be 
divided into simple four-sided zones and every three-dimensional region into six-sided 
zones. Grids can then be independently generated for each zone using existing grid- 
generation schemes. The zonal approach has the following advantages: 

1. The grid for any arbitrary region can be generated in a simple, straightfor- 
ward manner (since the zones are geometrically simple) . 
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2. Flow regions requiring grid refinement can be isolated in separate zones, 
and the required number of grid points can be introduced in these zones (if necessary 
this can be done adaptively as the solution progresses in time) . 

3. The zonal approach facilitates the use of different equation sets in differ- 
ent zones; hence, simpler equation sets can be used in certain regions of the flow. 
This may result in a saving of computer time. 

4. The zonal approach also facilitates a block-processing scheme wherein only 
the data corresponding to certain regions of the flow field are required to reside in 
the main memory of the computer; the remaining data can be stored on disk or tape. 
Theoretically, the block-processing technique permits the use of unlimited global grid 
sizes (computing speeds become a limiting factor in this case) . 


The division of a given region into zones introduces new boundaries in the cal- 
culation: the zonal boundaries. Since the grid for each region is generated inde- 

pendently, the grid lines of two adjoining regions may align (continuous grids) or 
may not align (discontinuous grids) with each other. Even in the case of continuous 
grids, a sudden change in grid spacing or grid-line orientation across the zonal 
boundary may give rise to discontinuities in the transformation metrics (metric- 
discontinuous grids). Figure 1 shows the different types of grids mentioned above. 


The following examples serve to illustrate some of the advantages of the zonal 
approach. Figure 2 shows the flow region associated with a combination of three air- 
foils. The region is multiply connected and hence difficult to discretize with a 
single grid system. The division of the flow field into 13 zones, as in figure 2, 
results in simple four-sided regions which can be discretized easily. Figure 3(a) 




b METRIC DISCONTINUOUS GRID 




c METRIC DISCONTINOUS GRID d DISCONTINUOUS GRID 

Fig. 1.- Types of grids used in finite-difference calculations. (a) Continuous, 
(b) Metric-discontinuous (abrupt change in grid line spacing) . (c) Metric- 

discontinuous (abrupt change in grid-line orientation). (d) Discontinuous. 
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shows a six-finned projectile, and figure 3(b) shows the grid used to discretize the 
flow region at the aft end of the projectile (ref. 1). Immediately after this axial 
station an entirely new type of grid (fig. 3(c)) is required because of the sudden 
termination of the body. This transition can be achieved by treating the plane corre- 
sponding to the aft end of the projectile as a two-dimensional (surface) zonal bound- 
ary separating two three-dimensional zones. It should be noted that the grid in 
figure 3(c) is itself composed of four two-dimensional grids separated by one- 
dimensional (line) zonal boundaries. The grids in figure 3(c) alleviate the geometric 
singularity problem associated with a simple polar grid and also facilitate the addi- 
tion of grid points to the central portion of the flow region (to accurately model the 
wake) . 

In order that information be transferred from one zone to another accurately, it 
is important to treat grid points on the zonal boundaries with care. The nonlinear 
nature of the Euler and Navier-Stokes equations permits solutions with discontinui- 
ties, such as shocks. It is imperative that the finite-difference scheme used for 
the calculation be conservative so that these discontinuities, when captured, assume 
the right strength and physical location. In a zonal calculation, it is important 
that the zonal boundaries also be treated in a conservative manner so that the dis- 
continuities can move freely across these boundaries. Results demonstrating the 
importance of maintaining conservation at zonal boundaries can be found in 
reference 2. 

A fully conservative zonal-boundary scheme, which permits the movement of dis- 
continuities across zonal boundaries with minimal distortion of the solution, is 
developed in reference 3. The scheme is designed for both discontinuous and metric- 
discontinuous grids and can be used in conjunction with first-order-accurate, 
explicit integration schemes. The scheme is stable, accurate, and is formulated in 
generalized coordinates. Results demonstrating these characteristics of the scheme 
are also presented in reference 3. The demonstration calculations include supersonic 
flow over a cylinder, blast-wave diffraction by a ramp, and the one-dimensional shock- 
tube problem solved on a two-dimensional grid. 

The first-order-accurate schemes used in reference 3 are insufficient to produce 
accurate results for a general class of problems. Hence, it is necessary to extend 
the zonal scheme of reference 3 so that it can be used with second-order-accurate 
integration schemes. The purpose of extending the explicit zonal scheme to work in 
conjunction with implicit schemes is twofold. The first reason is the increase in 
convergence rates and the consequent reduction in computing costs that can be obtained 
with the implicit schemes. Calculations with the Navier-Stokes equations require 
very fine grids in viscously dominated regions (to accurately resolve flow proper- 
ties) and relatively coarse grids in regions where viscous effects are negligible. 
Hence, mesh-cell areas in grids used for Navier-Stokes calculations can vary over 
several orders of magnitude. Navier-Stokes calculations with explicit schemes (with 
the attendant time-step restrictions) are impractical, and it is necessary to resort 
to implicit integration schemes. Since the ultimate objective of the zonal approach 
is to solve the unsteady Navier-Stokes equations in any given region, a first step in 
that direction is the implicit izat ion of the zonal-boundary condition. The advantage 
of using an implicit scheme, however, is not restricted to viscous calculations. The 
convergence rate can be expected to improve even in inviscid calculations where the 
mesh-cell size varies widely. The second reason for the implicitization of the zonal 
scheme is that it can then be incorporated into a number of existing analysis codes 
that use implicit integration algorithms. 
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The present study deals with the changes required to extend the zonal scheme of 
reference 3 so that it can be used with first-order-accurate and second-order- 
accurate implicit schemes, in particular with the Beam-Warming scheme (ref. 4) and an 
implicit form of the Osher scheme (ref. 5). The implicit zonal scheme and the inte- 
gration schemes used are described in the following sections. The zonal scheme is 
developed for the Euler equations cast in generalized coordinates. 

Results are presented for a cylinder in a supersonic flow. This calculation 
demonstrates the conservative property of the zonal scheme and also the considerable 
increase in convergence rates that can be achieved with an implicit integration 
scheme coupled with an implicit zonal scheme (relative to explicit integration and 
zonal schemes). A second calculation, consisting of a cylinder in subsonic flow, is 
also included. The calculation is performed on two grids, one rotating and the other 
stationary. This calculation demonstrates the time-accuracy of the scheme and also 
the feasibility of performing calculations on multiple zones of which some are moving 
relative to others. Zonal calculations in which zones move relative to each other 
will be frequently required in applications where some parts of the system move with 
respect to other parts, for example, the rotor-fuselage combination of a helicopter, 
turbines, and propeller-nacelle configurations. Finally, results are presented for 
the motion of a Lamb-type vortex through a zonal boundary. Pressure contours corre- 
sponding to the vortex remain distortion-free as they move through the zonal boundary. 
Since the exact position of the center of the vortex is known at all times (the prob- 
lem is formulated so that this information is known) , this calculation can be used to 
test the time-accuracy of the zonal scheme. The time-accuracy obtained in this case 
was found to be good. 

THE ZONAL-BOUNDARY SCHEME 


The explicit form of the zonal-boundary scheme is described in detail in refer- 
ence 3. In this section, the explicit zonal scheme is briefly outlined, and the 
necessary flux linearizations required to convert the scheme to an implicit scheme 
are then presented. 


Review of the Explicit Zonal Scheme 

Consider the unsteady Euler equations in two dimensions: 

CL + E + F = 0 (1) 

x t x y 

The vectors Q, E, and F are given by 


V 


pu 


pv 

pu 

E = 

p + pu 2 

F - 

puv 

pv 


puv 


P + pv 2 

e 


_(e + p)u_ 


_(e + p)v 


where p is the density; p is the pressure; u and v are the velocities in the x 
and y directions, respectively; and e is the total internal energy per unit volume: 
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Consider the two curvilinear grids used to discretize the flow region shown in 
figure 4. The line AB represents the zonal boundary that separates the two grids 
used to discretize the given region. Let £ and m be the indices used in the £ and 
n directions, respectively, in zone 1, and let j and k be the corresponding 
indices for zone 2. Note that the notation used in this study is different from that 
used in reference 3 (m and k increase in opposite directions in ref. 3). Let n 
represent the time-step for both zones. A superscript within parentheses will denote 
the zone to which a given quantity belongs; for example, Ar^ 1 ) denotes the marching- 
step size in zone 1. 



Figure 4.- Two-zone grid to illustrate zonal scheme in curvilinear coordinates. 

Establishing two independent-variable transformations (one for each grid) , 
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and applying these transformations to equation (1), the following is obtained: 


a (U 4 . a. pO — D i' — 19 

Q T d) E ? d) + F n U) ‘ 0 1-1,2 


( 5 ) 


where 
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( 6 ) 


= q/j^ 

= [C^Q + C (l) E + 5 (l) F]/J (l) 
t x y 

F (l) (Q,n (l) ) = [n^ l) Q + r/^E + i/ l) F]/J (l) 

l x y 

J<« - 5 ( «n<« - n (1) S (i) 
x y x y 

The notation E^^(Q»C^b and F^^(Q,n^^) is used to show the dependence of these 
quantities on the metrics of the transformation. Let the conservative difference 
schemes used to integrate equation (5) be given by 
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An 
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( 8 ) 


where E and F are numerical fluxes consistent with the transformed fluxes E and F. 
For an explicit integration scheme, these numerical fluxes are evaluated at the nth 
time-level, and for a fully implicit scheme they are evaluated at the (n + l)th 
time-level • 


The explicit zonal scheme consists of the following three steps: 

1. Integrate the dependent variables at grid points (of both the grids) that do 
not belong to the zonal boundary, using equations (7) and (8). 

2. Integrate the dependent variables at the zonal-boundary points of one of the 
zones (say zone 2 of fig. 4), using a scheme that conserves fluxes across the zonal 
boundary. 


3. Obtain the dependent variables at the zonal-boundary points of the other 
zone (say zone 1 of fig. 4) such that the dependent variables are continuous across 
the zonal boundary. 

The implementation of the first step of the zonal scheme is straightforward. 

The implementation of the second step is described below. Assume that the zonal- 
boundary points of zone 2 are to be updated using the finite-difference scheme of 
equation (8). This calculation requires the fluxes These fluxes have to be 
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calculated such that global conservation is maintained. A typical cell of a zonal- 
boundary point (j,l) is shown in figure 4, points RSTU. The points R and S are 
midpoints of the cells they lie in, and the points T and U are obtained as follows. 
The constant j lines of zone 2 are extrapolated into zone 1 to intersect the line 

CD (CD corresponds to m = mmax - 1/2 in zone 1 and to k = 1/2 in zone 2). The 

intersection points have the indices (j,l/2). Point T is midway between the points 
(j+1,1/2) and ( j ,1/2) » and point U is midway between points (j,l/2) and (j— 1,1/2). 

The global conservation property can be shown to be satisfied if the following rela- 

tionship is satisfied: 
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A close examination of equation (9) shows that each side of this equation is 
nothing but a discrete form of the line integral of the numerical flux F along the 
line CD in figure 4; the equation itself represents flux conservation across the 
zonal boundary. Equation (9) is only a necessary condition and is not sufficient to 
define the fluxes Fd '/ 2 a physically meaningful way. 


Assume that the 



are obtained by interpolating the 
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( 10 ) 


where the Nj^£ are interpolation coefficients and p and q define the set of points 
of zone 1 that will be used in the interpolation. We now describe a very simple way 
of obtaining the interpolation coefficients Nj £ such that equation ( 9 ) is auto- 
matically satisfied. Let the line CD in figure 5 correspond to the line CD in 
figure 4. The dots represent the grid points 



Figure 5.- Piecewise constant variation of the numerical flux F as a function of s. 
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of zone 1 and the crosses represent those of zone 2. A running parameter s is 
established along the line CD. The quantity s represents the distance of a point 
from the point C along the curve CD. Representative numerical values of 
^ mmax- 1/2 are pl° tte d on the positive y-axis. Assume a piecewise constant varia- 
tion of tne numerical fluxes Fji^mmax- 1/2 along CD, that is ^'sL,mmax—i /2 -*- s con- 
stant between s[i {/ 2 and s£+{/ 2 . Consider a point of zone 2, ( 3 , 1/2). Let U be 
midway between (j-1,1/2) and ( 3 , 1/2) and T be midway between ( 3 , 1/2) and (j+1,1/2). 

The F( 2 )/, are now calculated from 
J 9 * / 2 


p(2) CO _(2) 3 
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•j 

U 
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F^ , (Q , n n 1 ^ / ) ds 
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N F ^ 1 ^ (q n ( 1 ) \ 

j ,£ £,mmax-i/2 V 5 £ , mmax- 1 / 2 ' 


£=p 


(ID 


The end points in the interpolation p and q are shown in figure 5. The treatment 
of the points ( 1 , 1 ) and (jmax,l) and the calculation of the metrics used in the 
integration of the zonal points can be found in reference 3. 

The final step of the zonal scheme consists of updating the zonal-boundary points 
of zone 1 such that the dependent variables are continuous across the zonal boundary. 
This can be done very simply by an interpolation process. The dependent variables at 
the zonal-boundary points of zone 1 are obtained by interpolating the updated depen- 
dent variables at the zonal-boundary points of zone 2. The results in this study were 
obtained using a linear interpolation scheme. 


The Implicit Zonal Scheme 

The finite-difference equations (eqs . (7) and ( 8 )) can be made fully implicit by 
evaluating the numerical fluxes E and F at the (n + l)th time-level. However, this 
results in a system of nonlinear equations that in general need to be solved in an 
iterative manner. This problem of iteration (and the consequent increase in computing 
time) is usually circumvented by linearizing the numerical fluxes. The linearization 
process yields a system of linear equations that can be solved easily using a "block 11 
matrix solver. 

In one-dimensional problems, the block matrix resulting from the linearization is 
tridiagonal, and its inversion is computationally inexpensive. However, in two and 
three dimensions the matrix is still sparse, but it may have a large bandwidth 
(depending on the grid-point numbering system) . The inversion of the matrix in this 
case is computationally expensive. To overcome this problem, the left-hand side of 
the linearized equation is factored (ref. 4). Factorization results in multiple 
block-tridiagonal matrices on the left-hand side of the equation that need to be 
inverted sequentially. This approach has a disadvantage in that it introduces a fac- 
torization error that increases with increasing step size. In practice, the factor- 
ization error makes the scheme conditionally stable (in problems involving strong 
shocks that move rapidly through the grid, linearization errors could also lead to 
instabilities at large CFL numbers) . Factorization error also reduces the magnitude 
of the step size at which the fastest convergence is obtained. 
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Reference 5 presents a Newton-iterative scheme that retains the advantages of 
the linearization and factorization techniques. The iterative scheme for equation (8) 
can be written as 


I + xt [V C (A + ) P + A ? (A ) P ]J {i + [V n (B + ) P + A n (B _ ) P ]J[(Q (2) ) P+1 - (Q (2) ) P ] 


= -t(Q (2) ) P - (Q (2) ) n ] - At 


,(2) 


; ( 2 ) 




E v ” y _ E v_/ F v “ / _ p( 2 ) 

j+i / 2 ,k j- 1 / 2 ,k + j ,k+i / 2 j,k-i/2 


A? 


An 


( 12 ) 


where A + , A , B~*~, and B are Jacobian matrices given by 


^ = (A ± |A|)/2^ 
B 1 = (B ± |B|)/2 

> 

A = 3E/3Q 
B = 8F/3Q ^ 


(13) 


and A and V are first-order-accurate forward and backward difference operators. 

The variable p in equation (12) denotes the iteration number (Q p is an iterative 
update of Q n ) . When p = 0, = 0 n , and when equation (12) converges, Q p = Q n+1 . 

Note that the left-hand side of equation (12) becomes zero at convergence for every 
integration step. This implies that equation (8), with the fluxes evaluated at the 
(n + 1) th level, is satisfied at each integration step. A second point of interest is 
that the iterative scheme reverts to the noniterative scheme (of the type developed 
in ref. 4) when only one iteration is performed. 


In two dimensions, the iterative approach has the added advantage that the fac- 
torization error is driven to zero at each step if the iteration converges. It has 
also been found that the iterative approach enhances stability (ref. 5); that is, it 
is possible to use larger CFL numbers than with noniterative schemes. It should be 
noted that although the iterative approach permits the use of larger CFL numbers, 
there still exists a CFL limitation because at very large CFL numbers the iteration 
cycle may not converge owing to factorization error. This problem can be overcome by 
resorting to relaxation methods that do not require factorization (ref. 6). The 
primary disadvantage of the iterative approach is the added computation that the 
iterations require. This is offset to some extent by the larger step sizes that are 
possible with the iterative scheme. 


Noniterative factored schemes require grid lines to be continuous in order to 
maintain the block-tridiagonal structure of the implicit equations. Since grid lines 
are not continuous in zonal calculations, noniterative factored schemes are not the 
most efficient methods for solving the zonal, implicit equations. Iterative schemes 
are more adaptable to zonal calculations because approximations that restore the 
block-tridiagonal nature of the matrices can be made. Additionally, these approxima- 
tions can be corrected for during the iterative process. The details of such an 
approach are given below. 

In implicit zonal calculations it is required to linearize the numerical flux 
(Fj^'/ 2 ) n+1 (in addition to the interior fluxes). From equation (11) it is clear that 
this linearization will take the form 
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where , and , are Jacobian matrices given by 
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(15) 


The introduction of the AQj^mmax and A Q£^mmax-i into the system of linear equations 
of zone 2 destroys the block-tridiagonal nature of the system. In order to retain the 
block-tridiagonal nature of the system the following approximation is made: 


AQ, (1) . AQ< 2 > 

SI, umax 3,1 


£ = p 


(16) 


Equation (16) is a zeroth-order approximation which yields 
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(17) 


For a stationary grid, equation (17) reduces to 


AQ« - -1^ — AQ< 2 > 
£,mmax ,j(i) -) n J ’ 1 

' l,imax 


(18) 


Additionally, we make the approximation 


AQ^ = 0 

£ ,mmax-l 


(19) 


Equation (19) is an approximate Dirichlet-type of boundary condition. Equations (18) 
and (19) together decouple the calculations in the two zones (these approximations 
are corrected for in the iterative process). Substituting equations (18) and (19) 
into equation (14) yields 


<^:U >n+1 - < , i! > i/.> n + Z s T.l a l 

£=p 


mmax- 1/2 .-( 1 ) J 
' £,mmax' 


( 20 ) 
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For nonstationary grids, it is necessary to substitute equation (17) instead of equa- 
tion (18) into equation (14). The use of equation (20) instead of equation (14) to 
linearize the numerical fluxes at the zonal boundary restores the block-tridiagonal 
nature of the system. 


Equation (13) for the zonal-boundary point (j,l) of zone 2 now takes the form 


{i + 0- lv 5 (i + 


) F + A (A ) P ] 


'} 




ll+il . (D p + V M n+1 a p J- 1 ,,= (!). p+1 _ , B (2).p. 

) An tt Z-/ £ ,mmax“i/ 2 (1) n[ ^ ' (Q ) J 

£=p * £,mmax' J 


(J ) 

£,mmax y 


- RHS of equation (13) (21) 

Using the approximate linearization of reference 5, ot^nnnax-i/2 can now be written as 


a P = (B + ) P 

£,ramax-l/2 ' £,mmax' 


( 22 ) 


Note that if^the zone 1 boundary points were being integrated (instead of interpo 7 
lated) the B terms would come into effect. 


The interpolation coefficients 
can be replaced by for moving grids. This approximation is not required for 

stationary grids or for moving grids in which the movement is specified. The approxi- 
mation vanishes when the solution converges (or, if necessary, at each step, if the 
iterative process is permitted to converge at each step) . 


Clearly, equation (20) is only an approximation to equation (14). The inaccuracy 
in the linearization may lead to inaccuracies in the transient solution (the final 
converged solution, if one exists, will be the same). However, if the linearization 
given by equation (20) is used in conjunction with the iterative scheme described 
earlier (as in eq. (21)), time-accurate results can once again be obtained by iterat- 
ing equation (12) to convergence at each time-step. This is because the left-hand 
side is driven to zero at each time-step. The role of the approximate linearization 
given in equation (20), in the iterative approach, is to merely stabilize the calcula- 
tion during the iteration process (at CFL numbers greater than that allowed by an 
explicit scheme) . It should also be noted that equation (20) decouples the zone 1 
grid points from the zone 2 grid points (at the expense of multiple iterations). This 
decoupling permits a block-processing approach wherein only the data corresponding to 
a single zone, or even less, need to reside in the main memory of the computer. The 
information pertaining to the rest of the zones (except that used for the zonal- 
boundary condition) may reside on disk or tape memory. 

The procedure adopted in this study uses equation (20) in conjunction with the 
iterative scheme. In practice, only two iterations per time-step were required to 
perform stable zonal calculations up to CFL numbers of 50.0. At CFL numbers higher 
than 50.0, instabilities occurred at the surface boundary, and, hence, the zonal 
boundary condition could not be tested under such conditions. 

The implicit zonal scheme can be summarized in the following five steps: 


12 



1. Integrate the dependent variables at all the grid points of zone 2, using 
equation (12) in conjunction with the implicit zonal-boundary condition (eq. (21)) 
developed earlier (only one iteration) . 

2. Interpolate the newly obtained values of t(Qj^i) P+1 - (Qj^i) P l along the 
zonal boundary to yield a new set of values of [ (Q^nun ax ) p+1 ~ (Q^mmax^ P ^ • 

3. Integrate the dependent variables at the grid points of zone 1, using equa- 
tion (12) (only one iteration) and the most recent values of 

[ ^J^mmax )P+1 ” ^liax )P ^ (a Dirichlet-type Q f boundary condition). 

4. Interpolate the values of (Qj*i) P+1 to obtain the values of (Q£^ a x) P+1 
(discard the ones obtained as a result of the integration) . 

5. If the maximum value of the magnitudes of all [(Q^) P+1 - (Q^) P ] is less 
than the prescribed tolerance limit, go the next integration step; if not, go back to 
step 1 and reiterate* 

The scheme can be made time-accurate by choosing a sufficiently tight tolerance 
limit for the iteration cycle* For problems in which only the time-asymptotic solu- 
tion is required, step 5 may not be necessary. In the demonstration calculation 
included in this study (cylinder in a supersonic free stream) , step 5 was omitted 
and two iterations were performed at each time-step. 

Although a minimum of two iterations is required to perform stable calculations 
with the zonal technique described above, the second iteration need not be performed 
at all grid points. In the examples included in this study, it was sufficient to 
perform the second iteration only at the zonal-boundary points and at one row of 
points on each side of the zonal boundary. Thus, the extra expense incurred in using 
an iterative integration scheme as opposed to a noniterative integration scheme is 
typically less than 20% and decreases as the ratio of the number of zonal points to 
interior points decreases. A marginal increase in the number of steps to converge 
may be observed in some cases when the second iteration is performed only in the 
neighborhood of the zonal boundary. 


A Note on Conservation 

The iterative implicit zonal scheme (which uses eq. (20) at the zonal boundary) 
is conservative at each time-step only if the tolerance limit is sufficiently small; 
otherwise, the scheme is conservative only at convergence (if the problem has a time- 
asymptotic solution). However, if the zonal fluxes are correctly linearized, as they 
are in equation (14), then the scheme does not require iteration and is fully conser- 
vative at each time-step. The use of equation (14) would require the inversion of a 
sparse matrix that is not tridiagonal (it would also couple the dependent variables 
of the different zones) . The matrix could be either directly inverted or a relaxa- 
tion approach could be employed. If a relaxation scheme is used, then the 
approximate-factorization approach and its attendant problems can be discarded 
altogether. The relaxation approach to implicit zonal calculations will be presented 
in a separate article. 
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INTEGRATION SCHEMES 


The results of this study were obtained with the implicit, first-order-accurate 
and second-order-accurate Osher schemes and the Beam-Warming scheme. This section 
deals with the evaluation of the numerical fluxes for each of these schemes and with a 
slight modification that the zonal procedure requires in order to be applicable to the 
second-order schemes presented here. 


Numerical Fluxes for the Integration Schemes 


The numerical flux for the first-order-accurate Osher scheme is given by (ref. 7) 
1 


E j+i/2,k 2 + E(Q j+i,k ,C j+i/2,k ) 

- AE Xk’W’ W,k> + AE " (Q j,k’W’W,k )] 


(23) 


where 


AE (Q j ,k ,Q j+i ,k’^j+i/2 ,k } " 



dQ 


(24) 


A 

The numerical flux F^ can be obtained in a similar manner. 

The numerical flux for the second-order-accurate Osher scheme (ref. 7) can be 
summarized in the following expression: 

1 


W.k ' 2 + g( V».k’ 2 * * 5 j+l/2.k» 

-f [4E+ «J.k-Vi.k- 5 J-H/2.k> - “Vi.k^J.k’W.k” 


2 f AE ^j+i ,k’^j+2 >k’ ^j+ 1/2 ,k^ AE ^j,k ,Q j+i,k’ ? j+i/2,k^ ( ' 25 - ) 


where e - 1.0 for the Osher scheme. The Beam-Warming scheme (without smoothing) 

can be obtained by setting e = 0. For small values of e, we obtain the Beam- 
Warming scheme with smoothing, with the smoothing terms possessing an upwind flavor. 
Additional smoothing terms are not necessary for this form of the Beam-Warming scheme. 
The values of e that are usually used are 0.01 < e < 0*2. The upwind smoothing 

terms consist of both third- and fourth-order terms instead of the conventional purely 
fourth-order terms. Although the third-order nature of these terms does not affect 
the formal accuracy of the scheme, this fact should be kept in mind when viscous cal- 
culations are performed using the numerical flux of equation (25) , because for large 
values of e (say e > 0.1) the truncation error terms may become comparable to the 
natural viscous terms, thus reducing the accuracy of the calculation. However, small 
values of e (e < 0.05) are sufficient to stabilize the scheme, larger values being 
required only to damp out the oscillations near shocks. Since e can be varied over 
the region of interest without affecting the conservative property of the integration 
scheme, it is possible to use small values of e in viscously dominated regions and 
larger values of e in regions where shocks are captured. 
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The Beam-Warming scheme described above has the following advantages: 

1. The smoothing terms can be easily made conservative across the zonal boundary 
since they are part of the numerical flux. 

2. It is possible to transition to upwind schemes in a natural way in regions 
where upwind schemes are more appropriate than central-dif f erence schemes (regions 
where viscous effects are negligible). 

3. It provides a single framework in which both upwind and central-difference 
schemes can be represented. 

However, the disadvantage of using this form of the scheme is that it is compu- 
tationally more expensive than the conventional form (ref. 4). It should be noted 
that this additional expense is incurred only in viscously dominated regions where it 
is preferable to use the conventional scheme with fourth-order smoothing terms. 


Modifications to the Zonal Boundary Condition 

The procedure of calculating numerical fluxes to update the zonal boundary points 
requires a minor modification to be applicable to second-order accurate schemes that 
require a five-point difference molecule in each spatial direction. A close examina- 
tion of equation (25) will show that in order to calculate F^/^j we first need to 
determine the quantity AF + (Qj s „ ,Qj , x , n j } 3 / 2 ) • This is achieved’ simply by extrapo- 
lating the constant j lines of zone 2 (see fig. 4) to intersect the line 
m = inmax - 1 in zone 1 (or k = 0 in zone 2) instead of terminating the extrapola- 
tion at the line CD, which corresponds to k = 1/2 or m = mmax - 1/2. The values 
of Qj^o are interpolated from the values of Q^jnmax-i* A simple linear interpola- 
tion was used to obtain the results presented in this paper. A similar extrapolation 
(but this time into zone 2) is required to obtain the fluxes Fi^jnmax-i/a • The pro- 
cedure outlined above does not in any way affect the conservative property of the 
zonal-boundary condition. The additional flux terms that are required to make the 
integration scheme second-order accurate are once again perfectly balanced across 
the zonal boundary. 


RESULTS AND DISCUSSION 


Results are presented in this section for inviscid, supersonic flow past an 80° 
segment of a cylinder; for inviscid, subsonic flow past a full cylinder, and for vor- 
tex motion through a zonal boundary. The supersonic case demonstrates the increase 
in convergence rates that can be achieved with implicit zonal techniques and also the 
quality of solutions that the present zonal scheme yields. The subsonic case demon- 
strates the time-accuracy of the zonal scheme and the feasibility of performing 
calculations on zones that move relative to each other. The vortex calculation once 
again shows the time-accuracy of the zonal scheme and also the distortion-free motion 
of vortices through zonal boundaries. 


Cylinder in a Supersonic Free Stream 

The free-stream Mach number for this case of a cylinder in a supersonic free 
stream was chosen to be 2.0. The dependent variables at all grid points were 
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initialized to free-stream values. The finite-difference equations together with 
the various boundary conditions (including the implicit zonal-boundary condition 
developed in this study) were integrated to convergence. The boundary condition used 
at the surface of the cylinder is implicit and characteristic in nature. Details of 
its implementation are given in reference 8. 

Two— zone calculation— The grid used for the two— zone calculation is shown in 
figure 6. The discontinuity of the constant £ lines along the zonal boundary AB 
is clearly seen. The bow shock associated with this flow first appears at the surface 
of the cylinder and then moves outward through the zonal boundary to its converged 
position in zone 2. 


ZONE 2 


FREESTREAM 

BOUNDARY 



ZONAL BOUNDARY 
ZONE 1 


EXIT 

BOUNDARY 


CYLINDER 


SYMMETRY BOUNDARY 


Figure 6.- Grid for two-zone cylinder (supersonic) calculation. 
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Figure 7 shows the pressure contours 
obtained after 13 steps with the first-order- 
accurate, implicit Osher scheme. (The square 
symbols in figs. 7-10 represent the converged 
shock position predicted by another numerical 
approach in ref. 9.) Because of the large 
transients that occur during the first few 
steps of the calculation, the CFL number had to 
be gradually increased from 5.0 to 40.0. 

Table 1 gives the CFL numbers used at different 
stages of the calculation. The use of larger 
CFL numbers during the initial period of large 
transients resulted in the termination of the 
calculation, because of nonphysical solutions 
such as negative density. Figures 8 and 9 show 
pressure contours obtained after 19 and 24 inte- 
gration steps, respectively. The shock is seen 
to be passing through the zonal boundary in 
figure 8, and, in figure 9 the shock has passed 
through the zonal boundary. Figure 10 shows 
the pressure contours obtained at convergence. 
The contours in figure 10 are slope-continuous 
across the zonal boundary. This is because of 
the conservative nature of the zonal boundary 
condition and also because of the continuity of 
dependent variables that the zonal scheme 
enforces. 


Figure 11 shows the convergence history for 
the implicit zonal calculation (with the 
implicit, first-order-accurate Osher scheme and 
the implicit zonal-boundary condition) and for 
an explicit, zonal calculation (with the 
explicit, first-order-accurate Osher scheme and 
the explicit zonal-boundary condition) . The 
convergence criterion chosen was 


Ap < 5x10" 
1 1 max 


The implicit zonal scheme required 94 steps to 
converge, whereas the explicit zonal scheme 
required 2700 steps to converge; that is, the 
use of the implicit zonal scheme increased the 



A 


Figure 7.- Isobars after 13 inte- 
gration steps (first-order 
Osher scheme) . 


TABLE 1.- VARIATION OF CFL NUMBER 
WITH INTEGRATION STEP NUMBER 

Integration step (n) 

CFL number 

0 < n < 10 

5.00 

10 < n 

40.00 


overall convergence rate by a factor of 28.7. However, the implicit scheme with two 
iterations per time-step (the results shown in figs. 7—11 were obtained with two iter- 
ations per time-step) requires 4*8 times as much computing time per time-step as the 
explicit scheme. Hence, the actual computation cost was reduced by a factor of 6*0 
when the implicit zonal scheme was used. 


Although two iterations per time-step were used at all the grid points, only the 
zonal-boundary points and perhaps one or two adjacent rows require multiple iterations 
(two or more), because of the inaccurate linearization at the zonal boundary. A 
scheme that takes advantage of this should require less computation time per time- 
step than a scheme that iterates at all grid points. The two-zone calculation was 
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Figure 8.- Isobars after 19 integration 
steps (first-order Osher scheme). 


Figure 9.- Isobars after 24 integration 
steps (first-order Osher scheme) . 






-IMPLICIT OSHER SCHEME 
(IMPLICIT ZONAL BOUNDARY CONDITION) 


J I I I I 1 1 

0 400 800 1200 1600 2000 2400 2800 

NUMBER OF STEPS (n) 


Figure 10.- Isobars at conver 
gence (first-order Osher 
scheme . 


Figure 11.- Convergence history for the 
cylinder calculation (implicit and 
explicit first-order Osher scheme 
calculations) . 
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performed again, with two iterations being used only along the zonal boundary and the 
row adjacent to it. The CFL number was varied as in the previous case. The calcula- 
tion was found to be stable but the number of iterations for convergence increased to 
122 (a calculation with only one iteration at all grid points was found to be 
unstable). The computation cost (compared with that for the explicit scheme) was 
reduced by a factor of 7.3. Table 2 summarizes the computation times and convergence 
rates for the different approaches mentioned above. All computations were performed 
on a Cray-XMP. 


TABLE 2.- COMPUTATION TIMES FOR FIRST-ORDER-ACCURATE OSHER SCHEME 


Computation 

cost 

parameters 

Number of iterations 

per step 

2: Second at 2: Second on 

all points portion of grid 

Time per 100 steps, sec 

10.55 

20.18 

12.78 

Number of steps to 

Unstable 

94 

122 

converge 




Time to converge, sec 

Unstable 

18.97 

15.59 


The linearization of a fully implicit scheme results in the loss of some desir- 
able properties of the scheme (such as monotonicity in the case of the first-order- 
accurate Osher scheme). This can be seen clearly in figures 8 and 9. The oscilla- 
tions near the shock are seen even before the shock encounters the zonal boundary. 
These oscillations grow larger as the shock passes through the boundary and eventually 
vanish at convergence (fig. 10) as expected (ref. 5). The oscillatory behavior 
appears because only two iterations per time-step were used instead of the required 
number of iterations to converge the iteration process at each time-step. To demon- 
strate the restoration of the properties of the fully implicit scheme (when the iter- 
ations converge), the following test was performed. The solution obtained after 
22 integration steps (from the previous calculation) was used as the starting solu- 
tion and integrated through two steps. A CFL number of 40.0 with five iterations per 
step was used. The residual was reduced by approximately an order of magnitude at the 
fifth iteration. The pressure contours obtained as a result of this calculation are 
presented in figure 12. Clearly, the oscillations of figure 9 are absent in figure 12 
(figs. 9 and 12 can be directly compared since they show pressure contours at the 
same time-step) . The monotonicity property is restored because at convergence the 
difference equations corresponding to the fully implicit scheme are satisfied. 

Figure 13 displays the pressure contours obtained at convergence with the second- 
order-accurate Osher scheme. The contours are once again smooth and continuous across 
the zonal boundary. The second-order accuracy of the method results in oscillatory 
behavior of the contours near the shock, however, these oscillations are confined to 
a small region near the shock. The solution required 156 steps to converge when two 
iterations were performed at all the grid points and 155 steps when two iterations 
were performed only along the rows adjacent to the zonal boundary. The second-order 
scheme requires only about 25% more computation time per step than the first-order 
scheme. The total computation time to obtain the converged second-order solution is 
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Figure 12. Isobars after 24 integration Figure 13.- Isobars at convergence, using 

steps (the last two steps with 5 iter- the second-order Osher scheme, 

ations per step) , using the first- 
order Osher scheme. 

about 65/ more than that required by the first-order scheme. However, second-order 
accuracy is essential to building reliable, general-purpose codes. Results obtained 
with the first-order scheme have been presented merely to demonstrate the increase in 
convergence speeds and the decrease in computation costs (compared with the results 
of ref. 3) that are possible with implicit integration and zonal schemes. Table 3 
summarizes convergence speeds and computation costs for the second— order— accurate 
Osher scheme. 

Figure 14 shows the pressure contours obtained at convergence with the Beam- 
Warming scheme (e — 0.2). The oscillatory behavior of the solution near the shock is 
evident. Unlike the second-order Osher scheme, the oscillations are not confined to 
the Immediate vicinity of the shock, but extend as far as the zonal boundary. 

Although the oscillatory behavior seems to be associated with the zonal boundary, in 
reality it is not a result of the zonal boundary. This assertion is substantiated by 
the converged pressure contours shown in figure 15 obtained with the Beam-Warming 
scheme on the single— zone grid shown in figure 16. The oscillations appear once again 
and are similar to those in figure 14 (in fact the contours of figs. 14 and 15 are 
super imposable) • A summary of convergence speeds and computing times for the Beam- 
Warming scheme is given in table 4. 
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TABLE 3.- COMPUTATION TIMES FOR SECOND-ORDER-ACCURATE OSHER SCHEME 


Computation 

Number 

of iterations 

per step 

cost 

parameters 

i 2 

: Second at 

all points 

2: Second on 

portion of grid 

Time per 100 steps, sec 
Number of steps to 

12.96 

Unstable 

24.85 

156 

15.86 

155 

converge 

Time to converge, sec 

Unstable 

38.76 

24.58 



Figure 14.- Isobars at convergence, 
using the Beam-Warming scheme. 



Figure 15.— Isobars at convergence, using 
the Beam-Warming scheme in conjunction 
with the grid of figure 16. 
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Figure 16.— Grid for one-zone cylinder calculation with the Beam-Warming scheme. 


TABLE 4.- COMPUTATION TIMES FOR BEAM-WARMING SCHEME 


Computation 

cost 

parameters 


Number of iterations per step 


2: Second at 2: Second on 

all points portion of grid 


Time per 100 steps, sec 

12.96 

24.85 

15.86 

Number of steps to 

Unstable 

104 

149 

converge 




Time to converge, sec 

Unstable 

25.84 

23.63 
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Some of the oscillations shown in figures 14 and 15 can be removed, without 
affecting the accuracy of the solution, by adding additional, conventional fourth- 
order smoothing terms to the grid points in the vicinity of the shock. Such an 
attempt has not been made in this study. Also, the conventional Beam-Warming scheme 
requires less computer time per step than that shown in table 4. A conservative 
zonal scheme that uses the conventional Beam-Warming scheme will be discussed in a 
separate article. 

Four-zone calculation- A general-purpose, zonal Euler code should have the capa- 
bility of handling as many zones as necessary to cover the region of interest. In 
order to demonstrate the generality of the present zonal scheme and its applicability 
in a general-purpose, zonal Euler code, the region of interest for the cylinder was 
divided into four zones instead of into two zones as was done in the previous case. 
The zones and the grids for each zone are shown in figure 17. A wide variation in 
cell shapes and sizes can be seen across each of the three zonal boundaries. 

Figure 18 shows the pressure contours obtained at convergence with the implicit. 



Hill ill U HIHllll lZltI t 

t t 


ZONE 2 ZONE 1 

Figure 17.- Grid for four-zone cylinder 
(supersonic) calculation. 


Figure 18.- Isobars at convergence, 
using the first-order Osher scheme. 
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first-order-accurate Osher scheme. Small changes in contour line slopes can be seen 
at the zonal boundaries. This is because of the sudden change in accuracy caused by 
the abrupt transition in mesh clustering and will be less noticeable for higher-order- 
accurate schemes. Figures 19 and 20 show the pressure contours obtained with the 
second-order-accurate Osher and Beam-Warming schemes, respectively. Once again the 
contours are continuous across the zonal boundaries (and slope continuous unlike the 
contours in figure 18). The smooth transition of the shock from zone 2 to zone 4 
emphasizes the quality of solutions possible with the implicit zonal boundary scheme. 




Figure 19.- Isobars at convergence, using 
the second-order Osher scheme. 


Figure 20,- Isobars at convergence, using 
the Beam-Warming scheme. 


Cylinder in a Subsonic Free Stream (Two-Zone) 

The second type of calculation was performed for a complete cylinder in subsonic 
flow. The free-stream Mach number for this calculation was 0.35; figure 21 shows the 
two grids used in the calculation. All the grid points were initialized to their 
free-stream values and the governing equations together with the zonal boundary condi- 
tion were integrated to convergence (both the grids are stationary for this phase of 
the calculation) . The integration scheme used was the Beam-Warming scheme with 
e = 0.1. Second-order accuracy in time was achieved by using a three-point backward 
difference for the time-derivative on right-hand side of equation (12). The outer 
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ROTATING grid was then made to rotate at a constant 
OUTER angular speed. An integration scheme with- 

out any truncation error would result in the 
dependent variables being stationary in 
physical space and, consequently, contour 
levels of, for example, pressure and density, 
would also remain stationary (although in 
computational space, the dependent variables 
at a grid point of the outer grid would be 
constantly changing because of the changing 
physical location of the grid point) . 

Figure 22 shows the pressure contours 
obtained at each time-step, as the outer 
grid performed one rotation, plotted on the 
same figure. Two hundred and sixty inte- 
gration steps were used per rotation. Five 
iterations were performed at each step in 
order to reduce the magnitude of the resi- 
Figure 21.- Grid for two-zone cylinder due by approximately an order of magnitude. 

(subsonic) calculation. For an integration scheme of infinite accu- 

racy, the thickness of the bands would be 
zero. Since the integration scheme used is 
only second-order accurate in space and 
time, the thickness of the contour bands is 
finite. The thickness of the bands gives 
a qualitative measure of the time-accuracy 
of the integration scheme coupled with the 
zonal scheme. Further details of this cal- 
culation are given in reference 10. Refer- 
ence 10 also includes a demonstration cal- 
culation which shows that the finite 
thickness of the contour bands is a result 
of the accuracy of the integration scheme 
and not of any inadequacy in the zonal- 
boundary scheme. 

Vortex Motion through a Zonal Boundary 

This calculation consisted of a Lamb- 
type analytical vortex moving through a 
zonal boundary. It is possible to effect 
vortex motion either by superimposing a 
moving free-stream condition (in which case 
the vortex is convected along with the fluid 
at the free-stream velocity) or by keeping 
the vortex stationary and moving the grid. 
Figure 22.- Isobars at each step as the The two approaches yield identical results, 

outer grid performs one rotation. In this study, the vortex was initialized 

using the procedure given in reference 11, 

and then the grid was moved at a constant speed in a direction opposite to the direc- 
tion in which vortex motion was desired. Figure 23 shows the two-zone grid used for 
the calculation. Only the central portion of the flow field is presented in fig- 
ures 23-28, since the essential features of the vortex are contained in this region. 
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ZONE 1 ZONAL ZONE 2 

BOUNDARY 

Figure 23.- Grid for two-zone vortex calculation 



Figure 24.- Isobars for the vortex calculation (at initialization) 










The discontinuous nature of the grid lines at the zonal boundary is evident. The 
calculation was performed with the second-order-accurate Osher scheme and three iter- 
ations per step. Second-order accuracy in time was achieved as in the previous case. 

Figure 24 shows pressure contours at initialization. The solid core at the 
center of the constant-pressure circles in this figure and the following figures that 
show pressure contours is the analytical center of the vortex. Figure 25 shows pres- 
sure contours after 85 integration steps. The slope continuity of the contour lines 
across the zonal boundary is clearly seen. The center of the circles coincides with 
the analytical center of the vortex (this demonstrates the time-accuracy of the zonal 
scheme). Figure 26 presents pressure contours obtained after 180 integration steps. 
The vortex has moved entirely into zone 2. The contours are circular and are not 
distorted. However, because of asymmetric truncation errors, the center of the 
circles is slightly above the analytical center of the vortex. Figure 27 shows a 
continuous grid on which the calculation was performed once again. The results of 
this computation are shown in figure 28 (pressure contours after approximately 
140 integration steps). Once again a slight upward movement can be observed, thus 
demonstrating that this motion is not a result of an inaccurate zonal boundary 
condition. 


CONCLUSIONS 


An implicit, conservative zonal-boundary scheme that permits the use of discon- 
tinuous zonal grids is developed for the Euler equations. The scheme is an extension 
of an earlier zonal scheme that was designed for explicit integration schemes. The 
new implicit scheme can be used with first-order- or second-order-accurate integration 
schemes. The integration techniques used in this study are the implicit, first- 
order- and second-order-accurate Osher schemes, and the implicit Beam-Warming scheme. 
The special logic required to implement the zonal boundary conditions was found to be 
fairly simple to incorporate in existing codes. 

Results in the form of pressure contours are presented for inviscid supersonic 
flow over a cylinder with the associated bow shock. The calculation was performed on 
both two-zone and four-zone grid systems. The contours in both cases are observed to 
be continuous across zonal boundaries, and discontinuities were found to move freely 
and with minimal distortion through zonal boundaries. The zonal-boundary scheme was 
found to be stable even under severe test conditions such as strong discontinuities 
passing through the zonal boundaries at CFL numbers up to 40.0. At this CFL number, 
the use of the implicit zonal scheme increased the convergence rate by a factor of 28 
(over that of the explicit zonal scheme) , and the computation cost was reduced by a 
factor of 7. The results demonstrate the feasibility of efficiently using an implicit 
zonal approach with discontinuous grids in solving flow problems that involve discon- 
tinuous solutions . 

Results are also presented for subsonic flow past a cylinder calculated on two 
grids, one moving and the other stationary. These results demonstrate the time- 
accuracy of the zonal scheme and the feasibility of calculations on zones that move 
relative to each other. The ability to perform such calculations should prove useful 
in solving problems in which some parts of a system are moving relative to each other 
(e.g., rotor-fuselage combinations). A second problem that demonstrates the time- 
accuracy of the zonal scheme presented herein consists of a Lamb-type vortex moving 
through a discontinuous zonal boundary. The movement is distortion-free and 
time-accurate . 
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explicit schemes has been extended so that it can be used with existing 
second-order-accurate implicit integration schemes such as the Beam-Warming 
and Osher schemes. In the test case considered, the implicit schemes 
increased the rate of convergence considerably (by a factor of about 30 
over that of the explicit scheme) . Results demonstrating the time-accuracy 
of the zonal scheme and the feasibility of performing calculations on 
zones that move relative to each other are also presented. 
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